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Abstract 

We numerically investigate the fractal structure of two-dimensional quan- 
tum gravity coupled to matter central charge c for — 2 < c < 1. We re- 
formulate Q-state Potts model into the model which can be identified as a 
weighted percolation cluster model and can make continuous change of Q, 
which relates c, on the dynamically triangulated lattice. The c-dependence of 
the critical coupling is measured from the percolation probability and suscep- 
tibility The c-dependence of the string susceptibility of the quantum surface 
is evaluated and has very good agreement with the theoretical predictions. 
The c-dependence of the fractal dimension based on the finite size scaling hy- 
pothesis is measured and has excellent agreement with one of the theoretical 
predictions previously proposed except for the region near c ~ 1. 
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1 Introduction 



It could probably be one of the most annoying questions: "What is the observable 
of quantum gravity to check if the quantum theory of gravity is well defined ?" 
In two-dimensional quantum gravity the answer to this question is ready: "The 
fractal dimension of the quantum surface is the well-defined observable which can 
be analytically and numerically calculated." 

The importance of the fractal nature of two-dimensional quantum gravity was 
first recognized by KPZ [|TJ where the critical exponent is recognized to represent 
the fractal structure of the two-dimensional quantum surface. It has, however, been 
recognized later that more direct observable to detect the quantum fractal nature 
of space time would be the Housdorff dimension, or fractal dimension |fj. [J . 

The serious numerical study to measure the fractal nature of two-dimensional 
quantum surface was initiated by Agistein and Migdal ||]. They have proceeded the 
direct measurements of the fractal dimension of two-dimensional quantum surface by 
proposing a recursive sampling argorithm for c = model but could not observe the 
fractal nature of the quantum gravity in two dimension. The size of the triangulation 
was not large enough to observe the fractal nature of the quantum surface of c = 
model. The first numerical confirmation of the fractal structure of two-dimensional 
quantum gravity was carried out by Kawamoto, Kazakov, Saeki and Watabiki by 
using the recursive sampling argorithm for c = —2 model 0. In these numerical 
analyses the large size lattice triangulation (up to 5 million triangles) was necessary 
to confirm the fractal nature by the direct measurement of the fractal dimension. 
This numerical confirmation of the two-dimensional quantum space-time triggered 
wide varieties of numerical and analytic investigations of the fractal nature of two- 
dimensional quantum gravity. 

There are three analytic derivations of the c-dependence of Housdorff dimension 
or fractal dimension of two-dimensional quantum gravity coupled to matter central 
charge c, by Distler, Hlousek and Kawai |2|], Kawai and Ninomiya|| and later by 
Watabiki, Kawamoto and Saeki||. It was, however, pointed out that the measured 
fractal dimension of c = —2 model is very close to the third formulae given by 
Watabiki et. al. 0. In the meantime Kawai, Kawamoto, Mogami and Watabiki 
tried to understand the fractal nature of quantum gravity from the Matrix model 
point of view and succeeded to derive the transfer matrix of the quantum surface of 
two-dimensional pure gravity (c = 0)0. The formulation made it possible to derive 
the fractal dimension of pure gravity to be exactly (If = 4 which is consistent with 
the value of the first and third formulae. This analytic investigation of the random 
surface triggered further analytic and numerical investigations of two-dimensional 
gravityg i MM- 

Baby universe idea was proposed and proved to be useful to calculate string 
susceptibility numerically [|i~2"| , |i~3f . Then finite size scaling hypothesis was proposed 
later by being inspired by the analytic derivation of the two-point function of pure 
gravityjl4|]. The finite size scaling hypothesis made it possible to measure the fractal 
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dimension very accurately. Using this formulation we made systematic and the most 
reliable numerical measurement of the fractal dimension for c = —2 model fl5| , [16| . 
Then the very accurately measured fractal dimension of c = —2 model perfectly 
agreed with the theoretical value of the third formula, dp = 3.56±0.04(numerical) ~ 
3.561 • • • (theoretical). 

So far the numerical investigations of the fractal dimensions were carried out 
mainly for c = —2 and c = model and for several unitary series of conformal field 
theory |17| and several values in the region c > lRlS]. 



Here in this article we investigate the systematic investigation on the c-dependence 
of the fractal dimension of two-dimensional quantum surface. For c = —2, analytic 
formulae by the help of matrix model was available while for the other continuous 
value of matter central charge it was not obvious how to formulate the models to 
be convenient for the numerical study of the fractal dimension. It has, however, 
been known that the continuos central charge dependence can be accommodated 
by Q-state Potts model on the flat lattice. Here in this paper we reformulate the 
Q-state Potts model into the model which is a generalization of percolation cluster 
model, weighted percolation cluster model []i~9| 20| . We formulate the model on the 
dynamically triangulated lattice and extend to take non-integer value Q which is 
a function of c. By this model we can investigate c-dependence of the fractal na- 
ture of the two-dimensional quantum gravity coupled to matter central charge c. 
There have already been several calculations of critical indices of Ising model, three- 
state Potts models, and large Q values of Q-state Potts model coupled to quantum 



gravity 21, 22 



This paper is organized as follows: In section 2 we summarize the analytic deriva- 
tions of Hausdorff dimension and fractal dimension. In section 3 we explain the 
equivalence of Potts model and a weighted percolation cluster model which we call 
generalized Potts model. We then provide the definition of the generalized Potts 
model on the dynamically triangulated lattice. In section 4 we explain the details 
of the Metropolis argolism of the generalized Potts model. In section 5 numeri- 
cal results of the c-dependence of critical coupling constant, string susceptibility, 
and fractal dimension are shown. Conclusions and discussions are given in the last 
section. 



2 Analytic derivation of fractal dimension and pre- 
vious numerical results 

The fractal nature of the two-dimensional quantum gravity was first recognized by 
KPZ[|TJ. It was, however, not clear what kind of fractal it meant in the beginning. 
Serious analytic study on the fractal dimension of quantum gravity has been given by 
Distler, Hlousek and Kawai||, and Kawai and Ninomiya|3|] and later by Watabiki, 
Kawamoto and Saeki||. Here we summarize the analytic derivation of the fractal 
dimension. 
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In the derivation of the fractal dimension of two-dimensional quantum gravity 
coupled with matter central charge c, we use the formulation of Liouville theory in 



particular the formulation of conformal gauge given by DDK|^3|. We first summarize 
the main results of the formulation. 

Formally the continuum partition function for matter coupled to two dimensional 
gravity is given by 

Z(A) = Jvg 5^fd 2 x^-A^ Z M [g], (2.1) 

where Zu\g\ is a matter part of the partition function with gravitational background 
and A is the total area. 

Reguralized counterpart of the above partition function by dynamical triangula- 
tion is 

Zreg(A) = ^ Zm[G] S Na 2 A ~ Z M [Go], (2.2) 
G 

where N is the number of equilateral triangles and a 2 is the area of the triangle. 
G denotes a triangulation and Go is the typical triangulation which we select from 
the huge set of triangulations. The last approximate equality in eq.( [2.2|) is valid 
up to the normalization factor and if the selection of the typical surface is carried 
out by a correct procedure. Since the path integration of the metric is carried out 
after the selection of the typical surface, Go carries the information of the quantum 
fluctuation of space time effectively. 

Following David, Distler and Kawai (DDK) |23| , we obtain the gauge fixed version 
of the two-dimensional gravity coupled to matter central charge c with a conformal 
gauge = g^e*: 

Z(A) = J V^A fp [9]Zm[9] 6 (/ d 2 x fge a ^ - a) exp (^Ls^g]) , ( 2 .3) 
where App is the Fadeev Popov contribution. Sl [</>,<?] is the Liouville term given by 



S L [<b,g] = I d 2 x^g ( ~g^<f>d u <f> + <f>R) , (2.4) 



where we set the renormalized cosmological constant equal to zero for simplicity. a.\ 
appeared in eq. (|2.3|) is obtained from the following general formula: 



25- c - J(25 - c)(25 - c - 24n) 

«n = " ^ • (2-5) 

DDK have shown that the primary conformal field of weight n — 1 can be made 
Weyl invariant at the quantum level with a quantum correction: 

J d 2 x^e a "^ n , (2.6) 
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where <£>„ transforms as $ n \ge a= e^ n_1 ' ) ' J $ n |^. The term e" 1 * in the delta function 
of eq. ( |2.3D is needed to keep the world sheet volume to be A at the quantum level. 
Let us define an expectation value of an observable 0(g) by 

< 0(g) > A = Z(A)- 1 J V S (f)A F p[g]Z M [g]5 (J d 2 xfge a ^ - A 



' c — 25 
0(g,4>) exp ( 4g?r S L [<j),g] 



(2.7) 



2.1 Housdorff dimension and fractal dimension 

Here we define two types of critical exponents which specialize the fractal nature of 
the two dimensional random surface. 

Let us first define an intrinsic area A(r) of the random surface as a function of 
the mean square average size of the world sheet < r 2 > as viewed in the embedding 
space. We define Hausdorff dimension dn as 



A(r) = (V< r 2 >) dH . (2.8) 

It should be noted that this definition of Hausdorff dimension refers to the embedded 
space. 

Let us next define N(r) as the number of lattice points or number of triangles 
inside r steps from a marked site. A step on the original lattice is one link step of 
a triangle while a step on the dual lattice is one dual link (edge) step on the dual 
lattice. We define fractal dimension 4f as 

N(r)=r dF . (2.9) 

This definition of the fractal dimension characterizes the connectivity of the random 
surface, how the random surface is composed by the connection of triangles, and 
thus could be called connectivity dimension. 

Distler, Hlousek and Kawai|2| evaluated the mean square size of the quantum 
surface embedded in a D dimensional space by calculating the two-point Green's 
function of vertex operator 



1 2 

D <X >A 



dk 2 
1 / 



9 ln([ d'x.^gJxV) J d 2 x 2 ^g~(x7)e^ x ^- x ^ 

d 2 xi\]g(xi) j d 2 x 2 \lg(x 2 )(X(x 1 ) - X(x 2 )) 



k=0 



A 2 Z(A) \ 

A\ h: 



(A^oo), (2.10) 
where 7 S is the string susceptibility given by 



D-l - J (25 — D)(l — D) 
7, = ^ -• (2.H) 
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We can then obtain the first formula of Housdorff dimension as a function of D. 

= A> (2-12) 

I 7s I 

where D could later be identified as matter central charge c in two dimension. 

Let us next consider a derivation of fractal dimension using fermion as a test 
particle in the gravitational background following by Kawai and NinomiyaQ. The 
Lagrangian for matter fermion coupled to gravity can be given by 

i 

-yfgR + A^/g + eipi flip - meipip, (2.13) 



IQtt 

where A and e are, respectively, cosmological constant and the determinant of the 
vielbein in D-dimensions. In D = 2+e dimensions gravitational quantum corrections 
can be evaluated by the e-expansion formulation [Q. Then the fermion mass term is 
expected to acquire anomalous dimension via wave function renormalization of the 
matter fermion. 

Here we try to identify the anomalous dimension of the fermion mass term by the 
use of DDK formulation for Liouville theory. Under the scaling of the cosmological 
constant 

A ^/T D A, (2.14) 
the change in the Lagrangian can be absorbed by the field redefinition 

9„u -> /3V- (2-15) 

Since the vielbein transforms like square root of metric, the fermion kinetic term 
transforms as eipi pip — *> (3 D ~ 1 eipi pip. Then the scaling parameter can be absorbed 
by the field redefinition ip — > (3 <yl ~ D ^ 2 ip. The fermion mass term changes as meipip — > 
f3meipip, in particular ipip — > (3~ l ipip in D = 2. 

In two dimensions the fermion mass term m J cPxeipip is expected to have the 
form of eq. ( |2.6| ) after the introduction of the gravitational quantum correction. 
Since y/g and e acquire the same scale change, ipip can be identified as primary 
conformal field of $1/2 because of the same scale change: $1/2 \gp= /3 _1 $i/2 
and ipip —>■ (3~ l ipip. Then the Weyl invariant fermion mass term with the quantum 
correction is given by 

m [ d 2 xee ai ^ipip. (2.16) 



We now consider the following quantum average of the fermion mass term: 

mQ d 2 xee ai ^ipip^ , (2.17) 



where the quantum average < .. >a is defined in eq. ( |2.7|) . Under the constant shift 
of the conformal field 0^0 — 2/n/5/ai, the quantum average should be unchanged 
and yet the following relation holds: 

mQ d 2 xee a ^ipip^ = m(3~ 2 ^T Q d 2 x ee* 1 ^ ipip^j , (2.18) 
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where the following change of delta function is taken into account: 

5 (j d 2 Xy/g e ait/> - A\ -> f3 2 5 (j d 2 xy/ge a ^ - ft 2 A) . This relation suggests that 

/ n "l/2 \ 

the theory with two different sets of parameters (A, m) and I (3 A, m/3 Q i ) are 
equivalent. The two dimensional volume measured by the length scale of the fermion 



(2) 



field leads to another definition of fractal dimension d F 

,"1/2- " 



A ~ L 2 = [L ~j , (2-19) 
where 

dg>(c) = ^ = 2x v/ ^H + (2.20) 
with ct n given by the formula ( |2.5| ) . 



2.2 Fractal dimension from diffusion equation of random 
walk 

Here we provide yet another derivation of the fractal dimension by using the solution 
of diffusion equation on the random surface. 

We define the laplacian on the dynamically triangulated lattice. We first define 
adjacency matrix on the dynamically triangulated lattice. For a chosen typical 
surface Go we number the sites of the triangulated lattice. Then the (i, j) component 
of the adjacency matrix is defined as: = 1 if z-th site and j-th site are 
connected by a link, = if they are not connected by a link. It is interesting to 
note that (n, no) component of (K T )ij counts the number of possible random walks 
reaching from a marking site n$ to a site n after T steps. The Laplacian defined on 
the dynamically triangulated lattice is given by 

A L = 1 - S, Sy = —K^, (2.21) 

where qj is called coordination number and denotes a number of links connected 
to the site j. Sij is thus a probability of one step random walk from the site j to 
the neighboring site i. The diffusion equation on a triangulated surface Go with N 
triangles is given by 

d T ^°\T ]n ,n ) = A L (G Q )^°\T;n,n ), (2.22) 

where Ot is a difference operator in T and ^^ \T; n, n ) is a wave function of the 
diffusion equation and denotes the probability of finding the random walker at the 
site n after T steps from the starting site n . A solution of the diffusion equation 
can be obtained as ^^ \T; n, n ) = e TAi( ^ G ' ^(5 n) „, ), where (£„,„„) is iV-component 
vector with unit n entry. 
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We now consider the continuum limit of this diffusion equation. First of all we 
recover the lattice constant a. In taking continuum limit, the total physical area 
A = afNi is fixed and i — > oo; Oj — > 0, iVj — > oo is taken, where is the number 
of triangles and a 2 is the area of a triangle. In each step of the limiting process we 
select a typical surface Gi for the given number of triangles Ni, on which the lattice 
Laplacian Al(G{) of eq.( |2.21| ) is defined. Now the lattice version of the diffusion 
equation ( |2.22|) can be rewritten by 

\{^f'\T + a 2 ;x,x ) - Vf'\T;x,x )} = \a l {G^ *>(T; x, x )}, (2.23) 

a i a i 

where the location of the site x is measured with respect to the lattice constant a^. 
Thus we identify the dimension of T as that of area: dim[T] = dim[A\. In the con- 
tinuum limit the solution of the diffusion equation ( |2.23| ) is expected to approach 
the continuum wave function: (T; x, x ) — > \1/^ ?00 ' ) (T; x, x ). Numerically we 

approximate the limiting surface as the typical surface(Go) of the maximum size 
triangulation: G^ ~ G . As we have already noted in eq. ( |2.2j ), the metric integra- 
tion is effectively carried out for the equation ( |2.23j ) since we have chosen a typical 



surface. This means that the quantum effect is included for the wave function of 



eq.( 2.23 ). On the other hand the solution of the continuum counterpart of the diffu- 
sion equation: d T ^/(r; x, Xo) = A(g)^/(r;x,Xo) is still background metric dependent 
in general. Furthermore the dimensions of T and r may not necessarily be equal. 

Let us now define the comeback probability of random walk on the triangulated 
lattice and relate it with the continuum expression of Liouville theory as follows: 

G(T) = *f o) (T;* = x ) 

~ l^J d 2 x^^(r; x = x )J J Q d 2 x^j 

= d 2 x^e TA ^(0;x = x ) 



1 1 

A ~ T' 



(2.24) 



where < .. >a is the quantum average defined in eq.( p.7| ) and the last similarity 
relations are dimensional relations. We should remind of the fact that the metric 
integration is effectively carried out since we have chosen the typical surface Go for 
the wave function of the comeback probability. The initial wave function can be 
formally written as \l/(0;x = Xq) = lim x ^ xo 5 t (x — Xq)1/ \fg, where the regularized 
delta function is needed such as: 5 e (x — x ) = (l/V) x e/((x — x ) 2 + e 2 ). 

We next consider how to accommodate the Weyl invariance into the diffusion 
equation of random walk at the quantum level by using the formulation of Liouville 
theory. Let us consider the following quantity by Liouville theory: 



'J d 2 x v Tg^(T;x 
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Jd 2 x^/g^(0;x = x )j + r l^J d 2 x^fg A \P(0; x = x ] 



;2.25) 



where the solution of the diffusion equation is expanded by r. 

Taking a conformal gauge g^ v (x) = g^ v e^ x ' and introducing DDK arguments, we 

can rewrite the first and second terms of the eq.( |2.25| ) as 



J d 2 xy/g^{0;x = x )^ = /j d 2 x^Jg~ [-^5 e {x - x 



iX=XQ 

I A 



J d 2 x^gAV(0;x = x )^ = U d 2 xjg~e a -^ [X ±=8 e {x - x )} x=Xo 



A 

(2.26) 



where the term e" -1 ^ is introduced to keep the Weyl invariance of the second line 
of eq.( |2.26p in accordance with the arguments of eq. ( |2.6j ). Here it should be noted 
that A \I/(0; x = Xq) in the second term of eq. ( ]2.26| ) can be identified as a primary 
conformal field of $_i. 

Similar to the treatment for the Weyl invariance of the fermion mass term of 
eq. (|2.18j ), the quantum average of the comeback probability should be unchanged 
under the constant shift of the conformal field <fi — > <fi — 2lnf3/a\. In particular 
the second term of eq.( |2.25| ) should be unchanged and yet the following parameter 
change is generated: 

tQ d 2 x^e a - 1<t> A^(0;x = x )^ = t(3~ 2 ^T Q d 2 x ^fge a -^ A^(0; x = x )^ ^ 

(2.27) 

This relation suggests that the theory with two different sets of parameters (r, A) and 
(5 2 ~^~t, /3 2 A) are equivalent. Then we obtain the following dimensional relation: 



dimr = diim4 <*i . (2.28) 

We now point out that the expectation value of the mean squared geodesic 
distance is evaluated by the standard continuum treatment 

d 2 Xy/g {r(x, x )} 2 *(r; x, x ) 

d 2 x^/g {r(x,x )} 2 ^-^5(x - x ) + tA x -^5(x-x ) H ^ 

-At + 0(t 2 ), (2.29) 

which is related to the quantum version of the mean squared geodesic distance in 
the small r region 

< r 2 > 
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J2{r(x,x )} 2 ^°\T;x,x ) 

X 

J d 2 Xy/g J d 2 x y/g{r(x,xo)} 2 ^(r;x,xo)^ j \J d 2 x^fg 



A 



~ r ~ A Q i ~ T "i . (2.30) 

The last similarity relations are dimensional relations. Here we give the third defi- 
nition of fractal dimension 

j(3) 

A = (V< r 2 >) F (2.31) 

where 

4 3) (c) = -2^ = 2 x + ^fZ^ (2.32) 



with a n given by eq. 
2.3 Previous numerical results 

Serious numerical investigations of the fractal nature of two-dimensional quantum 
gravity (c = model) was initiated by Agistein and Migdal for c = modelQ. 
They proposed recursive sampling algorithm which partially use an analytic for- 
mula for tree diagrams and rainbow diagrams to generate planar Feynman diagrams 
corresponding to different types of two-dimensional sphere. They could't, how- 
ever, observe the fractal structure numerically for this model. The fractal nature 
of two-dimensional quantum gravity was numerically first confirmed by Kawamoto, 
Kazakov, Saeki and Watabiki by using the recursive sampling algorithm for c = —2 
model|J. It was recognized in this numerical analyses that large number of trian- 
gles is necessary to confirm the fractal nature by the parameterization of the fractal 
dimension given by fl2.9Q . In fact they needed 5 million triangles to confirm the 
fractal nature. It was then recognized that the number of triangles were not large 
enough to measure the fractal dimension of c = model. 

The numerical fractal dimension of c = —2 model was given in this numeri- 
cal analyses as dp ~ 3.55. After deriving the third formula of the fractal dimension 
df{c) of eq.flggg), we have recognized that the numerical value of the fractal dimen- 



sion is close to the theoretical value of the third formula d F (— 2) = (3 + vl7)/2 = 
3.561 

Except for the analytic derivation of the fractal dimension, one obtains the follow- 
ing analytic relations of comeback probability in (p.24[) and mean squared geodesic 
distance in ( |2.30| ) from the diffusion equation: 

1. (1) G(T)T ~ const., 

2. (2) < r 2 > = T 2 /^ 3) (~2) ~ T -o.56 ; 
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which can be numerically checked. These relations were numerically confirmed for 
c = —2 model in ||. 

An alternative analytic investigation of the fractal structure of the pure gravity 
(c = 0) were carried out by deriving transfer matrix of random surface by Kawai, 
Kawamoto, Mogami and Watabiki. They found out the beautiful scaling function 
p(L; D) which counts the number of boundaries whose boundary lengths lie between 
L and L + dL located at geodesic distance D measured from a marked point. It is 
evaluated by taking the continuum limit from the transfer matrix and disk amplitude 
of dynamical triangulation. The functional form of p(L; D) for c = model is given 
by " 

p(L; D)D 2 = -^(x- 5 / 2 + \x-^ + V)e-, (2.33) 

where x = L/D 2 is a scaling parameter. This quantity p(L; D)D 2 for c = model 
was measured numerically and had excellent agreement with the theoretical scaling 
function Q2.33|)Hll] . One of the important result of this analysis is that the fractal 
dimension of the c = model turns out to be dp = 4 which is consistent with the 
first and third formulae. 

Based on these theoretical and numerical evidences the third formula of the 
fractal dimension would be the correct formula for the c-dependence of the fractal 
dimension. In order to clear up the situation we started serious systematic and very 
accurate numerical study of c = — 2 model by using finite size scaling hypothesis fT5|, 
16fl . It was concluded that the measured fractal dimension from this analysis is 



dp = 3.56 ± 0.04 which is perfectly consistent with the theoretical value of the third 
formula. 



3 Generalized Potts model = Weighted percola- 
tion cluster model 

Potts model |H| was defined as a generalization of the Ising model [25]. Fortuin 



and Kasteleyn |20| showed that the Q-state Potts model is equivalent to a weighted 



percolation cluster model which we explain in this section. Their construction allows 
the Potts model to be generalized to nonintegral values of Q. We may call this model 
as generalized Potts model or equivalently weighted percolation cluster model. These 
models were originally formulated on the square lattice while we extend to formulate 
the model on a dynamically triangulated lattice. 

We define a planar y) 3 -graph Q dual to a triangulated lattice. Let a graph Qn 
have N vertices which are dual to triangles in the triangulated lattice. With each 
vertex i, we associate a spin that can take Q different values — 1, 2, • • • , Q. Two 
adjacent spins <7j and aj interact with interaction energy — J<5(<jj, <x,), where 
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Thus the Hamiltonian is 

H = -K^S(<T i ,<T j ), (3.2) 

where K = J/ksT can be reinterpreted as a coupling constant and the summation 
runs over all the pairs of nearest-neighbor vertices Then the partition function 

of this model is given by 



Z N {K;Q N )= ex P K ^2K a i^ a jj 



{a} on g N 



(3.3) 



where the cx-summation runs over all possible values {cr, = 1, 2, ■ • • , Q} for the spin 
variables <7i, a 2 , • • ■ , ctn on (?7v Thus there are Q N terms in the summation. 

It has been shown that the partition function (|3.3| ) can be expressed 

as a dichromatic polynomial ||27|| . In order to see this, let us expand the partition 
function as a product of terms associated with nearest-neighbor vertices. This can 
be worked out by using the relation [5(<Ti, crj)] m = 5(o"i, 0j) (m = 1, 2, ■ • •) as follows 

Z N (K; g N )=J2Ui 1 + v 6 ^ > ( 3 - 4 ) 

M (id) 

where we set v = exp(K) — 1. 

Let E be the number of the pairs of the nearest-neighbor vertices which we 
simply call edges on the graph Q^. It is equal to the number of original links in the 
triangulated lattice, i.e. 3N/2. Then the summand in eq. (|3.4j ) is a product of E 
factors. Each factor is the sum of two terms (1 and v 8{a iy <Tj)), so the product can 
be expanded as the sum of 2 E terms. Each of these 2 E terms can be associated with 
a bond-cluster-graph (from now on call this a cluster configuration) on Q^. Note 
that the term is a product of E factors, one for each edge. The factor for edge 
is either 1 or v 5{a^ crj): if it is the former, leave the edge empty, if the latter, place 
a bond on the edge with weight v 5(cjj, aj). Do this for all edges We then have 

a one-to-one correspondence between cluster configurations on Qn and terms in the 
expansion of the product in eq.( |3.4| ). 

Consider a typical cluster configuration, containing C connected components 
(1 < C < N), namely, clusters and b = Y^ =1 h bonds (0 < b < 3N/2). See fig.[l|. 
i-th cluster includes 6j bands and an isolated vertex may be regarded as a cluster 
with bi = 0. Then the corresponding term in the expansion contains a factor v b , 
and the effect of the delta functions is that all vertices within a component must 
have the same spin a. Summing over independent spins, it follows that terms gives 
a contribution Q°v b to the partition function Zjy. Summing over all such terms, i.e. 
over all cluster configurations, we have 

Z%[K-Q N )= £ Q c v b . (3.5) 

{cluster} 

on g N 

The summation is over all cluster configurations that can be drawn on Q^. The 
expression (|3.5|) is a dichromatic polynomial. Note that Q in eq.(p.5|) need not be 
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Figure 1: A fragment of a cluster configuration on a planar yj 3 -graph. An isolated 
vertex is regarded as a cluster. 

an integer. We can allow it to be any positive real number, and this can be a useful 
generalization of Q-state Potts model to non- integer real number Q. We may call 
this model as generalized Q-state Potts model or weighted percolation cluster model. 

Let us consider eq.( [3.5|) for a few particular values of Q. Firstly, we consider a 
model formulated on a given triangulated lattice. For the Q — > 1 limit, if 

we set v = p/(l — p) the partition function is 



V 



{cluster} 

y Ai-p) E - b - (3.6) 



i — p 



{cluster} 



Therefore, Z® 1 (K) becomes a sum over all possible bond percolation configurations 
with the correct weight where p is the probability of a bond being present . This 



result holds in any dimension and any lattice on which one defines the Potts models. 
Since the sum Y!, {duster} i n eq.( [3.6|) is the total probability and thus equal to 1, this 
model coupled to two-dimensional quantum gravity corresponds to the pure gravity 
model (c = 0). 

Next, let us examine for the Q — > limit. At the critical point of the Q-state 
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Potts model on the two-dimensional square lattice, it is known that v ~ Q 1 ! 2 |19| . 



In general, on a graph Q N we can assume v ~ Q a in the Q — > limit (0 < a < 1). 
Then the partition function Zjy °(.fr) becomes 

where / is the number of independent loops in a cluster configuration. We have used 
the Euler relation (See fig .[!]): 

b = N + l-C. (3.8) 



For < a < 1 the leading terms in eq. (|3.7| ) in the Q — > limit can be obtained by 
taking C — 1 (one cluster) and / = (no loops). These dominant configurations 
are just the spanning trees of the graph Qn |2(| p8| . Then each spanning tree 
configuration contributes to Z]y (K) with an equal weight. Therefore, this model 
coupled to two-dimensional quantum gravity is equivalent to the c = —2 scalar 
fermion model coupled to quantum gravity |J. 

Temperley and Lieb |29 used the result of Fortuin and Kasteleyn to prove the 



equivalence of the Potts model to the six-vertex or square-ice model, with staggered 
polarizations. Baxter, Kelland and Wu (BKW) |26l have later found a very elegant 
derivation of the result of Temperley and Lieb. They use a construction known 
as the BKW construction which makes many exact results obvious, including the 
critical temperature, self-duality and energy at criticality of the Potts model. 

Thus the Q-state Potts models are analytically solved and the relations with the 
conformal field theories are well known 0, |31]| . Q is related to central charge c in 
the following particular form: 



7T 



Q = 4cos 2 (— — ), c=l- " (3.9) 
\m + 1 / m[m + 1) 

The minimal unitary conformal field theories with central charge c between and 1 
correspond to integer m; m — 2, 3, 4, • • -. The generalization of Q to any positive real 
number corresponds to a continuous change of the central charge c, a generalization 
from minimal to non-minimal series of conformal field theories. 

Within the framework of dynamical triangulations the generalized Potts model, 
or equivalently the weighted percolation cluster models, coupled to two-dimensional 
quantum gravity is described by the following partition function: 

Z%(K)= £ J- £ Q c v b . (3.10) 

Gn^W 3 (T n )} {duster} 

where {<p 3 (T N )} denotes the set of <^ 3 (T/v)-graphs dual to triangulations T N of fixed 
topology (which we always assume to be sphere) and Sg N is a symmetry factor. 
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4 The Metropolis algorithm of the generalized 
Potts model 



We intend to evaluate numerically the fractal dimensions of two dimensional quan- 
tum gravity coupled to matter central charge c by the generalized Potts model, 
equivalently the weighted percolation cluster model formulated in the previous sec- 
tion ( p.lO| ). In the process of Metropolis updating we need double step updating: 
firstly we update a cluster configuration on a given graph Q N , secondly the graph 
Qn itself should be updated. 

Here we first formulate the Metropolis algorithm of cluster configuration on a 
given graph Q N . The equation (|3.5| ) leads us to choose the probability distribution 
of the generalized Potts model for a cluster configuration {C}k by 

P[{C}k] = 7*(K-C V (41) 

z N (K,y N ) 

where and b^ k ' are the number of cluster and the total number of edges = 
J2f=i b\ k ^ for the given cluster configuration {C}k, respectively. 

We need to define a transition function t[{C}k, {C}i\ for a transition {C}k — > 
{C}i, which satisfies ergodicity and the following detailed balance condition: 

P[{C} k ] t[{C} k , {CU = P[{C} Z ] t[{C} h {C} k ] . (4.2) 

Here we choose to use Glauber function |32| as the transition function 



t[{ c }k, {C}i] = j-^tt- , (4.3) 
1 + d Wik 



where 

t[{C} k , {C} t } _ mk = QSClkv 6 blk (44) 



t[{Ch,{C} k ] 

with SCi k = C w - C (fe) and 5b lk = 6® - 

For a given cluster configuration {C}k, our updating proceeds as follows: 

(1) We randomly pick up an edge on the graph Q N and change the edge by the 
following procedure and then the cluster configuration changes into {C}i. 

(2) We have to find out a change of the probability distribution 5Wik = Q sc "°v 5blk 
when we change the edge, where 5bik is the change in the total number of 
bonds. If the edge originally has a bond, remove the bond thus 5bik — — 1. If 
the edge originally doesn't have a bond, add a bond to the edge thus Sbik = +1. 
SC^ is the change in the number of clusters depending on the corresponding 
change of the edge. 
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(3) Next, we generate a pseudorandom number r uniformly distributed from to 
1 and change the edge following to the procedure (2) if and only if 



r < 



SWik 



(4.5) 



1 + SWi 



Ik 



This procedure ensures the transition {C}k — > {C}i with the correct proba- 
bility. We use the Glauber function for the transition function because of a 
faster convergence to the equilibrium distribution in our model. 

(4) Return to (1) unless the system is sufficiently equilibrated. 

By this updating, ergodicity and detailed balance condition can be ensured. 

Let us point out that it is nontrivial to evaluate 5Ca~ in the step (2). For a given 
cluster configuration {C}k, we pick up an edge on a graph Q N . Suppose the edge 
does not have a bond we add a bond on the edge with the probability of the step 
(3) and thus 5bik = +1 if the bond is added. For example we pick up the edge A-B 
which has vertices A and B in fig.^-(a) or the edge C-D which has vertices C and 
D in fig.0-(c), where those edges A-B and C-D does not have a bond. After the 
bond A-B and the bond C-D are added, (a) and (c) of fig.|| turn into (b) and (d) of 
fig.|2], respectively. In order to find 5Cn- we need to know if the both vertices of the 
edge belong to the same cluster or not before the bond is added. For example the 
vertices A and B in fig.|2|-(a) belong to the different clusters while the vertices C and 
D in fig.|2|-(c) belong to the same cluster. It is not time consuming to classify this 
difference numerically since we just need to know the data set of the collection of 
numbered vertices belonging to the same cluster. If both vertices originally belong 
to the different clusters then SCik = —1, while 5Cik = if they originally belong to 
the same cluster. It is thus numerically not difficult to identify 5Cik in the case of 
5b lk = +1. 

Let us suppose that we pick up an edge which already has a bond. Then we 
remove the bond with the probability of the step (3) and thus Sb^ = —1 if it is 
removed. We may pick up cluster configurations (b) and (d) of fig.|2] as particular 
examples which already has a bond at the edge A-B and C-D, respectively, and turn 
into (a) and (c) of figj^, respectively, after the bond A-B and C-D are removed. In 
order to find 5Cik we need to know the information if the both vertices of the edge 
belong to the same cluster or not after the bond is removed. The vertices A and B 
belong to the different cluster in (a) and thus SCik = +1, while the vertices C and D 
belong to the same cluster in (c) and thus 5Cik = 0. The crucial difference from the 
case Sbik = +1 to the case Sb^ = —1 is that the collection of the data set to classify 
the different cluster is not enough to differentiate if two vertices are in the same 
cluster or not after the bond is removed. The straightforward way to determine 
the connectedness of two vertices is to start from the first vertex and enumerate 
all vertices connected to it until either the second vertex is reached, or an entire 
connected region will be enumerated. This method is adequate if the clusters are 
small enough (tens of vertices), but if large clusters are involved (in the vicinity 
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(a) ' 






(d)^ 



Figure 2: Cluster configurations: (a),(b),(c),(d). Doded lines are edges without 
bond while solid lines are edges with bond. A and B are vertices in (a) and (b) 
while C and D are vertices in (c) and (d). Li denotes i-th surrounding loop. The 
arrows are pointers and compose the segments of the surrounding loops. 

of the critical point) the CPU time requirements can grow unreasonable. Since 
the connectedness, a non-local property, must be determined with every iteration 
to know 6Cik, it is important to find faster algorithm to evaluate 5Cik- In order 
to quickly determine the connectedness of large clusters, we have implemented an 
algorithm |33| using an auxiliary data structure based on the BKW construction 



26 



Let us reconsider a <p -graph Qn, with N vertices, together with its dual trian- 
gular lattice T N . If a bond is present on Q N , then its dual bond on T N is absent, and 
vice versa. The boundaries between clusters on and their dual clusters on T/v 
will form a collection of closed loops. Now we call this closed loops the surrounding 
loops. We then have the Euler relation, 



b + 2C = N + L. 



(4.6) 
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Figure 3: Three pointers surrounding each vertex. 

which relates the number of clusters and bonds to the number of surrounding loops 
L. 

By saving information of the surrounding loops as a data set, we transform the 
problem of determining connectedness of two vertices to the problem of determining 
whether two surrounding loop segments are part of the same surrounding loop or 
not. The surrounding loops are represented as a chain of pointers in the computer. 
A pointer is a memory location containing the address of the next pointer in the 
chain. There are three surrounding pointers represented by arrows for each vertex, 
as is shown in fig.[| We have shown surrounding loops composed of pointers for the 
figures (a), (b), (c), and (d) of fig.|[ For example in fig.||-(b) pointer Pi points to 
pointer P 2 , P2 points to P3, etc. In this manner the loops are represented by chains 
of pointers. Because of the (differential) Euler relation, Sbik + 25Cik = SLik, we can 
determine SCik if we can determine SLik, the change in the number of surrounding 
loops. 

Now, let us consider cases of removing a bond in the process of a Metropolis 
updating in Fig.|2|. The change; fig.|2|-(b) — > fig.Q-(a), illustrates the case in which 
removing the bond A-B will divide the surrounding loop Li into two surrounding 
loops Li and L 2 in figfj-(a), while the change; fig.0-(d) — > fig.^-(c), illustrates the 
case in which two surrounding loops L4 and Lq in fig.^-(d) will be joined into the 
loop L4 in fig|2|-(c) if the bond C-D is removed. In either case, two cuts must 
be made in chains of pointers and the four ends rejoined if the bond is removed. 
In the case of A-B bond in fig.^-(b), the P\ — > P 2 and the P4 — > P5 connections 
must be replaced by Pi — > P 5 and P 4 — > P 2 connections, respectively, in fig.Q-(a). 
Similarly in the case of the bond C-D in fig.Q-(d), the Pq — > P7 and the Ps — > Pg 
connections must be replaced by Pq — > Pg and Ps — > P7 connections, respectively, in 
fig.§-(c). By collecting the information of chains we can immediately conclude that 
SCik = +1 (for56^ = — 1) since Pi and P4 are in the same chain for fig.Q-(b) while 
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Figure 4: The flip move. This move is ergodic in the class of triangulations with a 
fixed number of triangles and fixed topology. 

Pq and P$ are in the different chains for figQ-(d) and thus 5Cik = (for 56^ = — 1). 
It is important to note that the imformation of the connectedness of the cluster 
configuration of (a) and (c) can be obtained by the data set of the chains of pointers 
of (b) and (d), respectively. It is numerically much easier to find if two pointers are 
in the same loop chain or not. 

In summary, in order to find 5Wik for the Metropolis algorithm of the generalized 
Potts model in Monte Carlo simulation, we first pick up an edge and find out 5bik 
depending on whether the edge has a bond or not. We first cut and rejoin our 
chains of pointers depending on addition or removal of the bond on the edge. Using 
the (differential) Euler relation, 5b[ k + 25Cik = SL [k , we evaluate 5C[k by judging 
whether two pointers attached to the different vertices of the given edge belong to 
the same chain or not. 

When we apply Monte Carlo simulations to quantum gravity coupled to the 
generalized Potts model or equivalently the weighted percolation clusters, we have 
to update the cluster configuration for a given triangulation and at the same time we 
have to update the triangulation for a given cluster configuration. The updating of 
cluster configurations on a given triangulation can be carried out just as described 
in the above. In order to update triangulations corresponding to the given (p 3 - 
graphs we use the standard flip-flop algorithm. We first choose two neighboring 
triangles randomly and flip the common link to generate a new triangulation. This 
flip move changes the triangulation locally as in fig.[|. This move is enough to make 
the process ergodic for the chosen class of triangulations; fixed number of triangles 
and fixed topology |J3|. 

In our simulations we avoid to generate configurations corresponding to tad- 
pole and self-energy graphs. In other words we consider the class of triangulations 
satisfying the following conditions: 
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Figure 5: Bond assignment after a flip move in a fragment of a cluster configuration. 
The bond is assigned to the new edge with 50 % probability. 

(1) No link has coinciding end. 

(2) Two sites may be connected by no more than one link (i.e. parallel links are 
forbidden) . 

In the updating of triangulations, if a triangulation in a forbidden class is generated 
we ignore the attempt to make the flip and instead choose a new link randomly and 
attempt a new flip. 

It is important to recognize that the flip procedure to generate random surface 
changes the cluster configuration. By the flip procedure the chosen edge which is the 
dual of the common link of the chosen neighboring triangles changes into a new edge 
which may or may not have a bond. In the present algorithm we generate a bond 
on the new edge with 50% probability. For example a cluster configuration in fig.[5] 
changes into two possible cluster configurations with an equal probability by the 
flip procedure. In this way a change in connectivity of two triangles concerned with 
updating entails a change in the assignments of neighbors and therefore a change 
in the weight function P[{C}fc] = Q c v b /Z®(K; Q N ) and then the transition 
function, 5Wj fc /(l + SWik), can be estimated according to the Metropolis scheme 
described in the above. 
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5 Numerical simulations and results 



The Monte Carlo simulations of generalized Potts models on dynamically triangu- 
lated surfaces are performed using the algorithm mentioned in the previous section. 
In our simulations we only generate planar y9 3 -graphs of spherical topology with no 
tadpole or self-energy configurations. And we perform simulations on graphs with 
the fixed number N vertices. 

We first define one sweep of our Monte Carlo simulation. In the process of 
updating for a given cluster configuration with a fixed triangulation one cluster- 
sweep means checking of the Metropolis procedure on each edge of the graph in a 
given order. On the other hand, in the process of updating dynamical triangulation, 
one triangulation-sweep means that we randomly pick up edges one after another 
with roughly 3iV/2 (total number of edges) times and try to flip them and at the 
same time we check to eliminate every triangulations containing tadpole or self- 
energy graph. 

After generating an initial cluster configuration, we thermalize our system by 
performing 2000 triangulation-sweeps and 1000 cluster-sweeps respectively. For the 
observables to be measured, this thermalization time is enough to thermalize our 
system according to the measurement of the autocorrelation time at criticality. To 
obtain a typical cluster sample configuration on a typical random surface background 
we perform 100 triangulation-sweeps and 20 cluster-sweeps. Then we measure sev- 
eral observables for the given configuration. We iterate these operations enough 
times until independent samples of a given number, which depends on the lattice 
size and on the observable to be measured, are obtained. 



5.1 The critical coupling constants 

It is a well accepted observation that lattice models with a second order phase transi- 
tion lead to corresponding continuum theories at the second order phase transition 
point. In particular minimal conformal lattice models lead to the corresponding 
continuum conformal field theory models at the critical point in two dimensions. 
It is well known that the Q-state Potts models on regular two-dimensional lattice 
correspond to the field theories of minimal unitary conformal series at the corre- 
sponding critical coupling constant K c . The correspondence is given in eq. (|3.9|) . It 
is analytically known that the Q-state Potts models make continuous (second order 
or even higher order) phase transition at the critical point for < Q < 4, while they 
have first order phase transition point for Q > 4. 

It is analytically not known if the nature of the phase transition may be changed 
when gravity coupled to the Potts models. For several examples of minimal unitary 
models, it is known that the order of phase transition can be raised to higher order 
when gravity is coupled |35, 36|. In the case of Q = 10 and 200, numerical results 



suggests a strong evidence in favor of continuous transitions |22|| . In our simulations 
we assume that the Q state Potts models coupled to gravity presented by eq. (|3.10| ) 
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have continuous phase transitions for < Q < 4 even at the non-unitary value of c. 

In order to locate the critical coupling by using finite-size scaling method, we 
investigate the percolation probability P(K) and the cluster size distribution n s (K) 
of percolation theory [B7J. Let us briefly summarize the physical meaning of the 
percolation probability P(K) and the cluster size distribution n s (K) by studying 
the simplest Potts model of Q = 1 case. The partition function of Q = 1 Potts 
model Z®~ l (K) is given by eq. ( |3.6|) . By the last expression of eq.( |3.6|) we can 
recognize that p with v = e K — 1 = p/(l — p) can be identified as the probability 
of a bond being present at an edge [pD| . When K — > 0,v — > and thus p — > 0, 
i.e. the probability of a bond being present at an edge is getting small, then the 
probability of finding a maximum cluster (the bonds of the maximum cluster reaches 
from one end to the other of the lattice extension) is expected to be zero. Let 
us define a quantity P^(K) = {the maximum cluster size}/iV, where the cluster 
size is defined as a number of bond of the cluster. Pn(K) is a probability of a 
vertex being on the maximum cluster, where N is the total number of vertices. 
P(K) = limjv^oo Pn(K) is called percolation probability. This quantity is the order 
parameter of the percolation transition and is expected to show the following critical 
singularity in the infinite lattice for \K — K c \ — > 

p[K ,[{K-K c f % K>K C 
K > \ 0, K < K c , { f 

where /3 P is the critical exponent associated to the percolation probability. 

In a finite lattice with the size N, the finite-size percolation probability Pn(K) 
can not vanish at any K > 0, and its behavior depends on the lattice size. In 
fig.p we show the size dependence of Pn(K) for Q = 2.5 and 0.6 as examples. 
As we can see, P^(K) with finite size dependence does not have sharp rise in 
contrast with eq. (|5.1|) but has milder rise with respect to K. In these simula- 
tions we have performed with the lattice sizes iV = 100, 200, 400, 800 and 1600 
for Q = 0.2, 0.4, 0.6, 0.8, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5 and 4.0 respectively. For each 
lattice size the number of independent samples is 5000. The range of the coupling 
constant K in which we measure the observables depends on the models, and they 
are shown in table [l]. So far as the behavior of Pn{K) is concerned, the order of 
phase transition is consistent with second order. 

For the behavior of Pn{K) near K — > K c we suppose the following scaling 
behavior based on the finite-size scaling hypothesis 



P N (K) = L^I V F P \{K - K c ) L 1 '"] , (5.2) 

where Q5.1p is assumed for the infinite system of iV — >• oo. L is the linear extension 
of the system and v is the critical exponent associated to the correlation length 
£,{K) ~ 1/\K — K c \ u . Since the total volume is proportional to N, we should make 
an identification iV = L dp with dp as fractal dimension. We may view £,(K) as 



an average geodesic distance. In order to extract K c using eq. (|5.2j ), we define the 
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K 



Figure 6: Size-dependence of Pn{K) for Q = 2.5 and 0.6 as examples. The sizes 
of the systems are N = 100 (the highest curve), 200, 400, 800 and 1600 (the lowest 
curve). 




Figure 7: The functions &n,n'(K), defined by eq. (p.3|) for all pairs of iVj's where 
N = 100, 200, • • • , 1600. The curves are Q = 2.5 and 0.6 as examples. 



following function $ NtN i(K) |3 



* m HP N {K)/P N ,{K)] 

] ^ JW] , (5.3) 

for a pair of size (N,N'). The functions 3>n,n'(K) and $>n',n"(K) for two different 
pairs of sizes (N, N') and (N', N") should thus intersect at K c , and the intersection 
point should yield —fipjvdp if corrections to finite-size scaling can be neglected. In 
fig.0 we plot the functions $jv,at'(-^) f° r au pairs of given sizes in the cases of Q = 2.5 
and 0.6. As we can see from fig.|6|, Pn{K) grows with the increase of Q, while the 
intersection point is relatively clear for larger Q. It is getting more difficult to find 
the intersection point for smaller Q. 

Let us next define a quantity n s (K) = p s (K)/s, where p s (K) is a probability 
of a vertex being on a cluster of size s. n s (K) can be recognized as a cluster size 
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Q 


range of K 


interval of K 


0.2 


0.25 ~ 0.82 


0.03 x 20 (points) 


0.4 


0.50 ~ 1.07 


0.03 x 20 (points) 


0.6 


0.70 ~ 1.27 


0.03 x 20 (points) 


0.8 


0.80 ~ 1.37 


0.03 x 20 (points) 


1.0 


0.83 ~ 1.40 


0.03 x 20 (points) 


1.5 


1.05 ~ 1.62 


0.03 x 20 (points) 


2.0 


1.23 ~ 1.65 


0.03 x 15 (points) 


2.5 


1.30 ~ 1.87 


0.03 x 20 (points) 


3.0 


1.42 ~ 1.84 


0.03 x 15 (points) 


3.5 


1.45 ~ 2.02 


0.03 x 20 (points) 


4.0 


1.55 ~ 1.97 


0.03 x 15 (points) 



Table 1: The ranges and intervals of the coupling constant K in which we measure 
the finite-size percolation probability P^{K) for various Q. 

distribution. This can be understood as follows: Suppose we have m s clusters with 
seize s, we can obtain a relation p s (K) = m s s/N and thus n s {K) = p s (K)/s = 
m s /N which is proportional to the cluster number and thus can be understood as a 
cluster size distribution. Using the cluster size distribution we define the so-called 
percolation susceptibility x(^0 

X (K) = J2's 2 n s (K), (5.4) 

s 

where the prime-sum means that the maximum cluster is omitted from the summa- 
tion. This quantity is the average number of bonds of a (finite) cluster and is ex- 
pected to show the following critical singularity in the infinite lattice for \K—K C \ — > 

(Tr . f (K-K e )-T, K>K C .... 
X{K) I {K c - K)^ , K < K c , (5 ' 5) 

where 7 P and 7^ are the critical exponents associated to the percolation susceptibility. 
Since these quantities P{K) and x(K) pl & y a major role in usual percolation theory, 
we use them as crucial quantities of determining the critical point of generalized 
Potts models. 

In a finite lattice with the size N, the finite-size percolation susceptibility xn{K) 
can not diverge at the critical point K c but reaches a maximum of finite height only. 
The magnitude of this maximum depends on the size of the lattice. In fig|| we show 
the size dependence of xn(K) for Q = 2.5 and 0.6 as examples. In these simulations 
we take lattice sizes N = 100, 200, 400, 800, 1600 and 3200 for various values of 
Q's. The range of the coupling constant K and the number of independent samples 
in which we measure the observables depend on the models and the lattice sizes, 
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1.4 1.45 1.5 1.55 1.6 1.65 1.7 0.6 0.7 0.8 0.9 1 1.1 

K K 



Figure 8: Size-dependence of xn(K) for Q = 2.5 and 0.6 as examples. The sizes of 
system are N = 100 (the lowest peak), 200, 400, 800, 1600 and 3200 (the highest 
peak). 



1 . 65 




Figure 9: The best linear fit to eq.( |5.7|) for Q = 2.5 and 0.6, where 5 = l/udp. 

and they are shown in tables |3|, f| and |5|. In a finite lattice there is so called pseudo- 
critical point K C (N) instead of the true critical point K c , where xn{K) reaches the 
maximum. The finite size scaling hypothesis for the correlation length |38| means 
that at the pseudo-critical point the correlation length coincides with the linear 
extension of the system, i.e. 

£(K e (N)) ~ L ~ N 1 ^ , (5.6) 

which, using £(K) ~ \K — K c \~ u , yields for iV — > oo 

\K C (N) - K c \ ~ N- 1/udF . (5.7) 

In order to extract K c we make three-parameter fit to eq. (|5.7| ). In fig.|^ we show the 
best linear fits for Q = 2.5 and 0.6. The error of K C (N) for each lattice size have 
been estimated using a polynomial fit to K C (N). Extrapolating to iV = oo yields 
the true critical coupling K c . 
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Figure 10: The best values of K c versus Q. The fat solid line is the fitting curve 
to the data points with polynomials of fifth-order in Q 1 ! 2 (gravity). We also show 
the theoretical critical coupling for the Q-state Potts models on the honeycomb and 
square lattice (flat). 



For large Q the critical couplings K c determined by using two methods as above 
are in agreement within the error. The discrepancy in the values is presumably due 
to the fact that larger lattices are needed before the methods converge completely. 
In fig.llO] we have plotted the best values of K c versus values of Q's. In this figure 



we have drawn the fitting function to data points with a polynomial of fifth-order in 
Q 1 / 2 . We also show the theoretical critical couplings for the Q-state Potts models 
on the honeycomb |4(| fy]] and square lattice (i.e. no longer coupled to gravity). 
As these two curves show a similar trend, this figure suggests the existence of exact 
solutions for the Q-state Potts models coupled to gravity. 

In the following subsections we use these values of the critical value K c obtained 



from the best-fit curve of fig. 10. 



5.2 The string susceptibility 

The string susceptibility exponent 7 S is one of the simplest quantities which charac- 
terize the fractal structure of quantum gravity. This quantity is introduced as the 
exponent of the subleading correction to the canonical partition function for random 
surfaces of fixed volume A 

Z{A) ~ e AcA A 7s " 3 , (5.8) 

for A — > oo, where A c denotes the critical cosmological constant. As we show 
in eq.(|2.2j), the total area A is proportional to the number of triangles N in the 
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Q 


K c 


0.00 


0.000 


0.05 


0.428 


0.20 


0.764 


0.50 


1.053 


0.80 


1.212 


1.00 


1.288 


1.50 


1.434 


2.00 


1.547 


2.50 


1.647 


3.00 


1.732 


3.50 


1.800 


4.00 


1.848 



Table 2: The critical coupling constants K c for various values of Q. 

regularized counterpart of the eq. (|5.8|) . In the case of pure gravity, Z(A = N) 
(a 2 = 1) is equal to the number of inequivalent triangulations with volume N. 

Physically the string susceptibility can be identified as an order parameter for the 
branching probability of random surface. This could be understood by the following 
relation: 

— !^ J\BBZ{B){A- B)Z(A- B) ~ A*', (5.9) 

where BZ(B) can be identified as a possible number of triangulation with a marked 
point on a triangulated surface of area B. Thus the left hand side of eq. ( |5.9| ) 
measures the average branching rate of the total surface branching into two parts p2| . 
For 7 S > the surface has tendency to branch more while for 7 S < the surface 
tends to be smooth. 

Using the analytical approach in a continuum framework, c-dependence of 7 S 
was first derived by KPZ Jl| and later rederived by DDK |23j by using conformal 
gauge formulation of Liouville theory, 



/ \ c-l-^/(25-c)(l-c) 
ls{c) = "—^ , (5.10) 

for two-dimensional quantum gravity coupled to the matter central charge c with 
spherical topology. 

In this subsection we investigate the string susceptibility exponent 7 S by mea- 
suring the distributions of so-called baby universes |12|, [13], [18| . It has already been 
pointed out that the numerical values of string susceptibility are in perfect agree- 
ment with the theoretical results ( |5.10| ) of Q-state Potts models for integral values 
of Q's. We then expect that the agreement will be perfect even for the non-integral 
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Figure 11: The measured string susceptibility 7 S versus central charge c and the 
theoretical curve. No logarithmic corrections are introduced in the fits. 



values of Q's. We intend to use the numerical investigations of 7 s (c) as the cross 
check of the critical values of K c calculated in the previous subsection. 

The branching probability of the surface with total area A into B and A — B 
is given by the integrand of eq. ( |5.9| ). The lattice counterpart of this branching 
probability is given by 



b N {B) 



3 (B + 1) Z(B + 1)(N-B + 1) Z(N -5 + 1) 



N 



B 

[!n -N 



~fs-2 



Z(N) 

[1«B< N), 



(5.11) 



where two baby universes are divided by single triangle in the current formulation 
of triangulations. We may call the smaller part of the minimum neck as a baby 
universe (minbu). In numerical simulations it is easy to find the shortest loops in a 
given triangulation and then to enumerate the area of the corresponding minbu's. 

The simulations were done with lattice sizes N = 1000 and 2000 for various 
values of Q. For each lattice size the number of independent samples is 100i^. The 
values of Q's and the critical couplings K C (Q) are shown in table ||. 

In order to extract 7 S we have fitted the distributions expressed in the form 



ln[b N {B)\ = ai + ( 7s - 2) In 



B 1 



B 

N 



+ 



a 2 
5 



(5.12) 



where a\ and a 2 are fitting parameters and the last term is a finite size correction 
term for small B [I3|. We have introduced a lower cut-off B c , because eq.( |5.11 ) is 
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Figure 12: The measured string susceptibility 7 S versus central charge c and the the- 
oretical curve. Logarithmic corrections are introduced in the fits for Q = 2.5, 3.0, 3.5 
and 4.0. 

only asymptotically correct, deviations can be expected for small B. Moreover we 
have cut large B part, because the distributions of minbu's are not universal for 
large B. 

The values j s (B c ) extracted from eq. ( |5.12|) approach exponentially to a limiting 
value for large B c . Thus in order to extract the limiting value 7 S we fit the values 
7 S (B c ) in such a form 



ls(B c ) 



Is - ol\ e 



-aiB c 



(5.13) 



In fig.[TI] we plotted the limiting value 7 S obtained by the above method versus 
central charge c together with the theoretical curve (|5.10| ). 

The reason why the results for c ~ 1 disagree with the theoretical curve is 
possibly due to the fact that logarithmic corrections are not yet introduced. It is 



well known |T2], |18[ that logarithmic corrections play an important role in the vicinity 
of c ~ 1. Thus we assume that the partition function has the following asymptotic 
behavior for large N 



Z{N) 



,A C JV 



N^~ 6 (lnAT) ( 



(5.14) 



where a is an additional parameter. The measured distributions b^{B) can now be 
fitted to the following parameterization: 



Oi + {is - 2) In 



B ( 1 



B 

N 



«2 



+ a ln[ln B ■ hx{N - B)\ + -± 

B 



(5.15) 
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for Q = 2.5, 3.0, 3.5 and 4.0. Then 7 S with the logarithmic corrections versus central 
charge c are plotted in fig.|T2"|. 

The string susceptibility 7 S with the logarithmic corrections for the various values 
of Q's are in very good agreement with the theoretical curve. We can then conclude 
that the values of critical coupling K C (Q) evaluated numerically in the previous 
subsection are correct within errors. 



5.3 The fractal dimension 

The most straightforward definition of the fractal dimension is given by 

N(r)=r djp , (5.16) 

where we count the number of triangles N(r) inside r steps from a marked triangle. 
In fact the fractal structure of the two-dimensional quantum gravity was first con- 
firmed in this way by Kawamoto, Kazakov, Saeki and Watabiki for c = —2 scalar 
fermion model [Q. In these analyses they needed 5 million triangles to obtain the 
reliable value of the fractal dimension. It was later recognized that finite size scaling 
hypothesis is very useful to evaluate the fractal dimension numerically and thus rela- 
tively small number of triangles is enough to obtain very accurate fractal dimension 



for c = —2 model Pql JIlGI. Here we use finite size scaling hypothesis to obtain the 
c-dependence of the fractal dimension. 

It has already been recognized numerically in [||] that the total perimeter length 
SAr(r) measured at geodesic distance r from a marked triangle grows 

s N ( r ) =r dp - 1 . (5.17) 

This fact triggered the investigation of analytic derivation of transfer matrix for two 
dimensional random surface of quantum gravity for c = model]?]]. It has then 
been recognized that two point function of random surface can be related to the 



measurement of s^ir) 0, [TJ], H3 



In the case of pure gravity (c = 0) "two-point function" with fixed volume A is 
defined by 

Sa(R) = / T>\g] 8( J d 2 x^g - A) i J d 2 x^ J d 2 x'^5(D g (x, x') - R) , 

(5.18) 

where D g (x, x') denotes the geodesic distance between the points labeled by x and x' 
measured with respect to g. Then Sa{R) dR is the average area of a spherical shell 
of thickness dR and radius R from a marked point on the manifold. We recognize 
that the lattice triangulation version of Sa(R) corresponds to Sjy(r). According to 
the numerical result, we expect to have a relation 

S A (R) ~ R dF ~ l for R ~ 0. (5.19) 
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Figure 13: The best linear fit to eq. (|5.26| ) for Q = 2.5 and 0.5 as examples, by the 
decay of the peak of s^(r)/N with N. 

In case of pure gravity (c = 0) Sa{R) was calculated analytically 

S A {R) = R 3 f{R/A*), (5.20) 



where f(x) ~ e 2:4/3 for large x [Q, 13]. This analytic result is consistent with the 



calculation of the fractal dimension dp = 4 for c = model derived by transfer 
matrix formalism 0. This result, however, strongly suggests the following scaling 
hypothesis for general matter central charge c coupled two dimensional quantum 
gravity: 

S A (R) = A l - l / dF F(R/A l ' dF ) , (5.21) 

where 

F(x)~x dF -\ x<l. (5.22) 

This beautiful scaling behavior was observed with a high accuracy in the very sys- 
tematic numerical analyses of the fractal structure of two-dimensional quantum 



gravity coupled to c = — 2 matter [TJ|, jig] . Here we assume that the scaling hypoth- 
esis works in general for the general matter central charge —2 < c < 1. 

In our numerical simulations we define a geodesic distance as a "triangle dis- 
tance" instead of link distance on a triangulation for saving the CPU time. The 
triangle distance denotes the shortest path along neighboring triangles between the 
two separate triangles. Therefore the triangle distance is equal to the "edge dis- 
tance" in the dual <y9 3 -graph. Then the discretized analogue of S A (R), sjy(r) with 
A = Na 2 (a 2 = 1), consists of all triangles with triangle distance r measured from a 
marked triangle. Corresponding to the above continuum description we expect the 
following behavior for Stv(^): 

sn(t) ~ iV 1-1 ^ F(x) , x = j^, (5.23) 
and we expect F(x) to behave as x dp ~ l for small x. 
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Figure 14: The best linear fit to eq. (|5.27| ) for Q = 2.5 and 0.5 as examples, by the 
decay of the inverse average radius l/(r)jv with N. 

In order to measure correlation function SAr(r) the simulations were carried 
out with lattice sizes ranging; N = 100, 200, 400, 800, 1600, 3200 and 6400 for 
Q = 0.0, 0.05, 0.2, 0.5 and 0.8. For Q = 1.0, 1.5, 2.0, 2.5, 3.0, 3.5 and 4.0 we have 
performed the simulations with lattice sizes ranging; N = 100 ~ 6400 and 12800. 
For each lattice size ranging; N = 100 ~ 1600 the number of independent samples 
are lOOi^, and we choose 20 initial random vertices (rv) on each configuration. Then 
the number of independent samples is 50K x 40(rv) for N = 3200, 20i^ x 100 (rv) 
for N = 6400 and 10K x 200 (rv) for N = 12800, respectively. The values of Q's 
and its critical couplings K C (Q) are shown in table |2]. 

In order to extract the fractal dimensions dp from the scaling relation ( |5.23| ), It is 



crucial to introduce the so-called shift parameter a ||T4| , |15| , |I~6|1 which accommodate 
the finite size effects. This parameter a is considered as the leading order correction 
in the scaling variable x = r/N 1 ^ dp . This is reasonable from the point of view that 
the shortest distances may include lattice artifacts, and thus we can not expect 
an exact agreement with the continuum formulae. On the other hand, in order 
to incorporate the higher order corrections in the scaling variable x we need to 
introduce a second shift parameter b [ TA\ : N 1 / d " F — > jV 1//rfF + b. In this way we are 
led to a "phenomenological" scaling variable x 



x 



r + a 



N^dF + h ■ 



(5.24) 



In the present simulations with small lattice sizes (at most iV ~ 10 3 ) we need to 
introduce the parameter b, however, it is numerically very hard to determine the 
optimal values of three parameters a, b and dp in the case of the gravity coupled 
to matter. Although such optimal values of parameters are obtained, they keep the 
dependence on lattice size. 

We have extracted the fractal dimensions (Jp by the following method. First we 
use the decay of the peak of SN{r)/N with iV (this is similar to the use of finite-size 
scaling in the case of the percolation susceptibility). From eq.( |5.23| ) we expect the 
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following scaling behavior: 

s N (r)/N ~ N- l/dF . (5.25) 

We have divided iV's into three successive divisions such as iV's ((100,200,400), 
(200,400,800), etc.), and we have performed to fit s N (r)/N = aiN~ 1/dF for each 
successive divisions by two parameters. In this way the dp dependence on N becomes 
clear. Then we further assume the following scaling behavior: 

\l/d F (N) - l/d F (oo)| ~ N~ a , (5.26) 

where we can extrapolate 1/^(00) by a linear fit to eq.(5.26). This is shown in 



fig. 13 for Q = 2.5 and 0.5 as examples, where the estimation of error is based on a 



non-linear fits to Sm(r)/N. The values of dp's obtained in this way for the various 



values of Q's are plotted in fig.fTE Three theoretical curves given by the formulae; 



(p. 12 ), (|2.20|) and (|2.32| ), are shown to be compared. It is clear that the formula 



.321) is closer to the numerical values of the fractal dimension. 



Secondly, in order to make use of the whole information of S7v(r) we have used 
the average radius (r)jv of universes with volume N 

= T7 E r s N(r) ~ N 1/dF . (5.27) 

^ V r=0 

We can then expect l/(r)jv ~ N~ l l dp , where d F can be determined in the same way 
as above. In fig.[l4] we show the corresponding linear fits to eq.( |5.27| ). The values of 



dp's obtained in this way for the various values of Q's are plotted in fig.|TB 
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These two independent results are consistent with each other and support the 
theoretical prediction d F '(c). 



6 Conclusion and Discussions 

In this article we have shown the results of numerical investigations of the fractal 
dimension of two-dimensional quantum gravity coupled to the matter central charge 
c for —2 < c < 1. The c-dependence of the matter central charge is introduced 
by reformulating Q-state Potts model into the model which can be identified as a 
generalization of percolation cluster model, weighted percolation cluster model, on 
the random lattice. In this formulation Q can be generalized into non-integer value 
and thus continuous c-dependence is realized and then we have called this model 
simply as the generalized Potts model. Since the model has a percolation cluster 
feature, we have formulated a new Metropolis algorithm to generate clusters on 
dynamically triangulated surface. 

The c-dependence of the critical coupling K c is not known theoretically. We 
have evaluated the c-dependence of K c by measuring percolation probability and 
percolation susceptibility with the help of scaling hypothesis. The c-dependence of 
the critical coupling has similar behaviour as those of flat lattice. It is then very 
natural to expect that there is a theoretical solution of theoretical coupling K c in two- 
dimensional quantum gravity coupled to matter central charge c. The order of the 
phase transition is assumed to be second or even higher order at the critical point and 
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we have not observed any evidence against this assumption. The string susceptibility 
is measured by the baby universe technique and has excellent agreement with the 
theoretical curve. It is recognized that the next leading logarithmic correction is 
important to improve the agreement with the theoretical prediction in particular for 
the region near c ~ 1. Several parameters were needed to incorporate the finite size 
effects. The string susceptibility measurement is carried out as a cross check of the 
critical values of K c and the result was perfectly consistent with the independent 
measurement of the critical exponent in the above. 

The c-dependence of the fractal dimension is measured based on the two-point 
function of quantum surface with the finite size scaling hypothesis. We needed two 
extra parameters to accommodate the finite size corrections to the scaling parameter. 
Measurements are carried out by two methods; the decay behaviour of the peak of 
two-point function and the average radius of universes. The results agree well each 
other. The c-dependence of the fractal dimension is excellent agreement with the 
following theoretical prediction except for the region c ~ 1: 



4 3 )(c) = 2 x V^+ V^C 
V25 - C + y/1 - c 

We consider that the deviation of the agreement of the numerical values of the 
fractal dimensions near c ~ 1 from the theoretical values is possibly due to the fact 
that the size of the lattice is not large enough to observe the possible discrete jump 
from the fractal dimension c < 1 of eq. ( |6.1| ) to the value of branched polymer phase 



It is interesting to measure the change of fractal dimension very accurately in 
this delicate region near c ~ 1. The theoretical curve of the c-dependence of the 
fractal dimension has infinite slope at c = 1 and then turns into imaginary value. We 
conjecture that the fractal phase of two dimensional quantum gravity c < 1 turns 
into branched polymer phase in 1 < c. We expect that there is a discrete jump of 
the fractal dimension at c = 1. In order to measure this discrete jump numerically 
we may need huge number of triangles in the simulation. 

It is interesting to compare with the measurement of the string susceptibility 
in three-dimensional simplicial gravity. In three dimensions tetrahedron is the fun- 
damental simplex which corresponds to the triangle of two-dimensional quantum 
gravity. In three-dimensional dynamical triangulation the gravitational constant 
can be a free parameter and plays a role of central charge c of two-dimensional 
quantum gravity. It was measured that the string susceptibility changes from neg- 
ative region to the positive region where branched polymer phase is expected |47| . 



Here again the phase change from the fractal phase to the branched polymer phase 
is expected. For realistic higher dimensional simplicial quantum gravity it would 
be important to understand the phase change such as the fractal-branched polymer 
phase change. 
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Table 3: The lattice sizes, ranges of the coupling constant K and the number of 
samples in which we measure the finite-size percolation susceptibility xn{K) for 
Q = 0.2, 0.4, 0.6 and 0.8. 
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20000 




800 


1.05 


r^j 


1.32 


0.03 


X 


10 


(points) 


20000 




1600 


1.13 


r-^j 


1.31 


0.02 


X 


10 


(points) 


20000 




3200 


1.14 


r-^j 


1.34 


0.02 


X 


11 


(points) 


10000 


1.5 


100 


1.01 




1.39 


0.02 


X 


20 


(points) 


20000 




200 


1.13 


r^> 


1.41 


0.02 


X 


15 


(points) 


20000 




400 


1.19 


r^> 


1.47 


0.02 


X 


15 


(points) 


20000 




800 


1.28 


rsj 


1.46 


0.02 


X 


10 


(points) 


20000 




1600 


1.31 


rvj 


1.49 


0.02 


X 


10 


(points) 


20000 




3200 


1.33 


rsj 


1.51 


0.02 


X 


10 


(points) 


10000 


2.0 


100 


1.10 


CNJ 


1.60 


0.02 


X 


26 


(points) 


20000 




200 


1.20 


CNJ 


1.60 


0.02 


X 


21 


(points) 


20000 




400 


1.28 


CNJ 


1.64 


0.02 


X 


19 


(points) 


20000 




800 


1.37 


C*NJ 


1.61 


0.02 


X 


13 


(points) 


20000 




1600 


1.41 


/*NJ 


1.63 


0.02 


X 


12 


(points) 


20000 




3200 


1.45 




1.63 


0.02 


X 


10 


(points) 


10000 


2.5 


100 


1.26 


P*J 


1.64 


0.02 


X 


20 


(points) 


20000 




200 


1.37 




1.65 


0.02 


X 


15 


(points) 


20000 




400 


1.42 




1.70 


0.02 


X 


15 


(points) 


20000 




800 


1.50 


rvj 


1.68 


0.02 


X 


10 


(points) 


20000 




1600 


1.53 




1.71 


0.02 


X 


10 


(points) 


20000 




3200 


1.54 




1.72 


0.02 


X 


10 


(points) 


10000 



Table 4: The same of table | for Q = 1.0, 1.5, 2.0 and 2.5. 
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Q 


# triangles 


rang 


e of K 


interval of K 


# samples 


3.0 


100 


1.35 




1.73 


0.02 


X 


20 


(points) 


20000 




200 


1.46 


r*j 


1.74 


0.02 


X 


15 


(points) 


20000 




400 


1.50 


r*j 


1.78 


0.02 


X 


15 


(points) 


20000 




800 


1.59 


r*j 


1.77 


0.02 


X 


10 


(points) 


20000 




1600 


1.60 


r*j 


1.80 


0.02 


X 


11 


(points) 


20000 




3200 


1.62 


r*j 


1.80 


0.02 


X 


10 


(points) 


10000 


3.5 


100 


1.40 


r*j 


1.78 


0.02 


X 


20 


(points) 


20000 




200 


1.53 


r-^j 


1.81 


0.02 


X 


15 


(points) 


20000 




400 


1.58 


r-^j 


1.86 


0.02 


X 


15 


(points) 


20000 




800 


1.66 


r-^j 


1.84 


0.02 


X 


10 


(points) 


20000 




1600 


1.68 


p*J 


1.86 


0.02 


X 


10 


(points) 


20000 




3200 


1.69 


r-^j 


1.87 


0.02 


X 


10 


(points) 


10000 


4.0 


100 


1.42 




1.90 


0.02 


X 


25 


(points) 


20000 




200 


1.56 




1.90 


0.02 


X 


18 


(points) 


20000 




400 


1.65 




1.93 


0.02 


X 


15 


(points) 


20000 




800 


1.73 


r*j 


1.91 


0.02 


X 


10 


(points) 


20000 




1600 


1.74 


r*j 


1.92 


0.02 


X 


10 


(points) 


20000 




3200 


1.75 


r*j 


1.93 


0.02 


X 


10 


(points) 


10000 



Table 5: The same of table ^| for Q = 3.0, 3.5 and 4.0. 
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